#include <fstream>
#include <iostream>
#include <iomanip>
using namespace std;

#include "plant_2DSOS_RDMD.h"
// --------------------------------------------------
plant_2DSOS_RDMD::plant_2DSOS_RDMD(double idt,int iL,int iMode, double iE_D,double iE_S,double iE_N)
{
	// Options
	LatticeSize = iL;
	mode = iMode;
	assert(2*iMode<=iL);
	E_D = iE_D;
	E_S = iE_S;
	E_N = iE_N;
	FLAG_CAL_MODE = 1;

	// 2 Inputs: #0 = T, #1 = W
	int input_num = 2;
	int input_dim[2] = {1,1};
	// 4 outputs:
	//  Output #0~3: surface profile
	//  Output #4: h
	//  Output #5: r2
	//  Output #6: m2
	int output_num = 7;
	int output_dim[7];
	output_dim[0] = (mode+1)*(mode+1);
	output_dim[1] = (mode+1)*(mode+1);
	output_dim[2] = (mode+1)*(mode+1);
	output_dim[3] = (mode+1)*(mode+1);
	output_dim[4] = 1;
	output_dim[5] = 1;
	output_dim[6] = 1;
	block::initialize(input_num,input_dim,output_num,output_dim,idt);
	
	cov_z1 = outputs[0];
	cov_z2 = outputs[1];
	cov_z3 = outputs[2];
	cov_z4 = outputs[3];

	// States
	lattice = new int*[LatticeSize];
	lattice_stat =  new int*[LatticeSize];
	for(int i=0;i<LatticeSize;i++)
	{
		lattice[i] = new int[LatticeSize];
		lattice_stat[i] = new int[LatticeSize];
	}
	p_migr = new double[LatticeSize];
	p_desorp = new double[LatticeSize];

	ef_cos_table = new double*[mode+1];
	ef_sin_table = new double*[mode+1];
	int m,i;
	double SQRT_TWO = sqrt(2.0);
	double SQRT_PI = sqrt(pi);

	for(m=0;m<=mode;m++){
		ef_cos_table[m] = new double[LatticeSize];
		ef_sin_table[m] = new double[LatticeSize];
		for(i=0;i<LatticeSize;i++){
			if(m == 0){	
				ef_cos_table[m][i] = 1.0/SQRT_PI;
			}else{ 
				ef_cos_table[m][i] = SQRT_TWO/SQRT_PI*cos(2*m*pi/LatticeSize*i);
			}
			ef_sin_table[m][i] = SQRT_TWO/SQRT_PI*sin(2*m*pi/LatticeSize*i);
		}
	}
	reset();
}
// --------------------------------------------------
plant_2DSOS_RDMD::~plant_2DSOS_RDMD()
{
	int i;
	for(i=0;i<LatticeSize;i++)
	{
		delete []lattice[i];
		delete []lattice_stat[i];
	}
	delete []lattice;
	delete []lattice_stat;
	delete []p_migr;
	delete []p_desorp;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::reset()
{
	block::time = 0.0;

	// Inputs
	T = 0;
	W = 0;

	// States
	int x,y;
	for(y=0;y<LatticeSize;y++){
		for(x=0;x<LatticeSize;x++){
			lattice[y][x] = 0;
		}
	}
	for(y=0;y<LatticeSize;y++){
		for(x=0;x<LatticeSize;x++){
			init_local_stat(x,y);
		}
	}

	// Outputs
	int m,n,index;
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D_mode(m,n);
			cov_z1[index] = 0.0;
			cov_z2[index] = 0.0;
			cov_z3[index] = 0.0;
			cov_z4[index] = 0.0;
		}
	}
	h  = 0.0;
	r2 = 0.0;
	m2 = 0.0;
}
// ---------------------------------------------------
void plant_2DSOS_RDMD::init_local_stat(int x,int y)
{
	lattice_stat[y][x]=(lattice[y][(x+LatticeSize-1)%LatticeSize]>=lattice[y][x])
		+(lattice[y][(x+LatticeSize+1)%LatticeSize]>=lattice[y][x])
		+(lattice[(y+LatticeSize-1)%LatticeSize][x]>=lattice[y][x])
		+(lattice[(y+LatticeSize+1)%LatticeSize][x]>=lattice[y][x]);
	return;
}
// ---------------------------------------------------
void plant_2DSOS_RDMD::init_migr_rate()
{
	int i;
	double k_m;

	k_m=2*T*KB2H;

	for(i=0;i<4;i++)
	{
		MIGR_RT[i]=k_m*exp(-(E_S+i*E_N)*EV2KB/(T+0.0));
	}
    MIGR_RT[4]=0;
	return;
}
// ---------------------------------------------------
void plant_2DSOS_RDMD::init_desorp_rate()
{
	int i;
	double k_m;

	k_m=2*T*KB2H;

	for(i=0;i<5;i++)
	{
		DESORP_RT[i]=k_m*exp(-(E_D+i*E_N)*EV2KB/(T+0.0));
		//DESORP_RT[i]=0;
	}
	return;
}
// ---------------------------------------------------
void plant_2DSOS_RDMD::init_p_migr()
{
	int x,y;

	migr_ttl=0.0;
	for(y=0;y<LatticeSize;y++)
	{
		p_migr[y]=0.0;
		for(x=0;x<LatticeSize;x++)
		{
			p_migr[y]+=MIGR_RT[lattice_stat[y][x]];
			
		}
		migr_ttl+=p_migr[y];
	}
	return;
}
// ---------------------------------------------------
void plant_2DSOS_RDMD::init_p_desorp()
{
	int x,y;

	desorp_ttl=0.0;
	for(y=0;y<LatticeSize;y++)
	{
		p_desorp[y]=0.0;
		for(x=0;x<LatticeSize;x++)
		{
			p_desorp[y]+=DESORP_RT[lattice_stat[y][x]];
			
		}
		desorp_ttl+=p_desorp[y];
	}
	return;
}
// ---------------------------------------------------
void plant_2DSOS_RDMD::update(double sysTime,double iT,double iW)
{
	if(time>sysTime)
		return;
	
	double rdm;
	while(time<=sysTime)
	{
		if(iT!=T)
			set_temperature(iT);
		if(iW!=W)
			set_ads_rate(iW);
		
		rate_ttl = ads_ttl+migr_ttl+desorp_ttl;
		rdm = RandRand()*rate_ttl;
		if(rdm<ads_ttl)
			adsorption();
		else if(rdm<=migr_ttl+ads_ttl)
			migration();
		else
			desorption();

		time += OneMCTime(rate_ttl);
	}
	if(FLAG_CAL_MODE){
		h2mode();
	}
	outputs[4][0] = getH();
	outputs[5][0] = getR2();
	outputs[6][0] = getM2();
	return;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::update(double sysTime)
{
	update(sysTime,inputs[0][0],inputs[1][0]);
}
// --------------------------------------------------
void plant_2DSOS_RDMD::getSnapshot(const string FileName)
{
	ofstream fout(FileName.c_str());
	fout << "LatticeSize  = " << LatticeSize << "\n";
	fout << "LatticeSize  = " << LatticeSize << "\n";
	fout << "Time = " << time << "\n";
	int x,y;
	for(y=0;y<LatticeSize;y++)
	{
		for(x=0;x<LatticeSize;x++)
		{
			fout << setw(5) << lattice[y][x];
		}
		fout << "\n";
	}
	fout.close();
}
// --------------------------------------------------
double plant_2DSOS_RDMD::getM2()
{
	int x,y;
	double sum = 0;
	double m2;
	for(y=0;y<LatticeSize;y++)
	{
		for(x=0;x<LatticeSize-1;x++)
		{
			sum+=(lattice[y][x]-lattice[y][x+1])*(lattice[y][x]-lattice[y][x+1]);
		}
		sum+=(lattice[y][LatticeSize-1]-lattice[y][0])*(lattice[y][LatticeSize-1]-lattice[y][0]);
	}
	m2 = sum/LatticeSize/LatticeSize;
	return m2;
}
// --------------------------------------------------
double plant_2DSOS_RDMD::getH()
{
	int x,y;
	double h_avg,sum = 0;
	for(y=0;y<LatticeSize;y++)
	{
		for(x=0;x<LatticeSize;x++)
		{
			sum+=lattice[y][x];
		}
	}
	h_avg = sum/LatticeSize/LatticeSize;
	return h_avg;
}
// --------------------------------------------------
double plant_2DSOS_RDMD::getR2()
{
	int x,y;
	double R2;
	double sum=0;
	double h_avg;
	
	h_avg = getH();
	for(y=0;y<LatticeSize;y++)
	{
		for(x=0;x<LatticeSize;x++)
		{
			sum+=(lattice[y][x]-h_avg)*(lattice[y][x]-h_avg);
		}
	}
	R2 = sum/LatticeSize/LatticeSize;
	return R2;
}
// --------------------------------------------------
double plant_2DSOS_RDMD::getTime()
{
	return time;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::set_temperature(double T_set)
{
	T = T_set;
	init_migr_rate();
	init_desorp_rate();
	init_p_migr();
	init_p_desorp();
	k_m2d = exp(-(E_D-E_S)*EV2KB/(T+0.0));
	return;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::set_ads_rate(double W_set)
{
	W=W_set;
	ads_ttl=W*LatticeSize*LatticeSize;
	return;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::add_remove_atom(int x, int y, int num)
{
	int xpos[5],ypos[5];
	double site_rate_mig, site_rate_desorp;
	int i;

	GetNeighbors(x,y,xpos,ypos);
	for(i=0;i<5;i++)
	{
		site_rate_mig = MIGR_RT[lattice_stat[ypos[i]][xpos[i]]];
		site_rate_desorp = DESORP_RT[lattice_stat[ypos[i]][xpos[i]]];
		migr_ttl  -= site_rate_mig;
		desorp_ttl -= site_rate_desorp;
		p_migr[ypos[i]] -= site_rate_mig;
		p_desorp[ypos[i]] -= site_rate_desorp;
	}

	if((lattice[y][x]==0) && (num==-1))
	{}
	else
	{
		lattice[y][x]+=num;
	}

	for(i=0;i<5;i++) 
	{
		init_local_stat(xpos[i],ypos[i]);
	}

	for(i=0;i<5;i++)
	{
		site_rate_mig=MIGR_RT[lattice_stat[ypos[i]][xpos[i]]];
		site_rate_desorp=DESORP_RT[lattice_stat[ypos[i]][xpos[i]]];
		migr_ttl+=site_rate_mig;
		desorp_ttl+=site_rate_desorp;
		p_migr[ypos[i]]+=site_rate_mig;
		p_desorp[ypos[i]]+=site_rate_desorp;
	}
	return;
}
// --------------------------------------------------	
void plant_2DSOS_RDMD::GetNeighbors(int x, int y, int *xpos, int *ypos)
{
	xpos[0]=x;
	xpos[1]=x-1;
	xpos[2]=x;
	xpos[3]=x+1;
	xpos[4]=x;
	if(xpos[1]<0) xpos[1]+=LatticeSize;
	if(xpos[3]>=LatticeSize) xpos[3]-=LatticeSize;

	ypos[0]=y;
	ypos[1]=y;
	ypos[2]=y-1;
	ypos[3]=y;
	ypos[4]=y+1;
	if(ypos[2]<0) ypos[2]+=LatticeSize;
	if(ypos[4]>=LatticeSize) ypos[4]-=LatticeSize;

	return;
}
// --------------------------------------------------	
void plant_2DSOS_RDMD::adsorption()
{
	add_remove_atom((int)(RandN1()*LatticeSize),(int)(RandN1()*LatticeSize),1);

	return;
}
// --------------------------------------------------	
void plant_2DSOS_RDMD::migration()
{
	int x,y;
	int xpos[5],ypos[5];
	int i,n;
	int st[4][2];

	migr_site(&x,&y);
	if(lattice_stat[y][x]==4)
		return;

	GetNeighbors(x,y,xpos,ypos);
	
	// Count the number of neighbors with lower height
	for(i=1,n=0;i<5;i++)
	{
		if((lattice[y][x]-lattice[ypos[i]][xpos[i]])>0)
		{
			st[n][0]=ypos[i];
			st[n][1]=xpos[i];
			n++;
		}
	}

	if(n==0) return;
	
	n = (int)(RandN1()*n);
	add_remove_atom(x,y,-1);
	add_remove_atom(st[n][1],st[n][0],+1);

	return;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::desorption()
{
	int x,y;
	migr_site(&x,&y);
	add_remove_atom(x,y,-1);
	return;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::migr_site(int *migr_x, int *migr_y)
{
	int y,x;
	double rdm;

	rdm = RandN1()*migr_ttl;
	for(y=0;rdm>=0;y++)
		rdm-=p_migr[y];
	y--;
	rdm+=p_migr[y];
	for(x=0;rdm>=0;x++)
		rdm -= MIGR_RT[lattice_stat[y][x]];
	x--;

	*migr_x = x;
	*migr_y = y;
	
	/* FOR DEBUG
	int event[5];
	for(i=0;i<5;i++)
	{
		event[i] = 0;
	}
	if(lattice_stat[y][x]==0)
		event[0]++;
	else if(lattice_stat[y][x] == 1)
		event[1]++;
	else if(lattice_stat[y][x] == 2)
		event[2]++;
	else if(lattice_stat[y][x] == 3)
		event[3]++;
	else if(lattice_stat[y][x] == 4)
		event[4]++;
	*/
	return;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::desorp_site(int *desorp_x, int *desorp_y)
{
	int y,x;
	double rdm;

	rdm = RandN1()*migr_ttl;
	for(y=0;rdm>=0;y++)
		rdm-=p_desorp[y];
	y--;
	rdm+=p_migr[y];
	for(x=0;rdm>=0;x++)
		rdm -= DESORP_RT[lattice_stat[y][x]];
	x--;

	*desorp_x = x;
	*desorp_y = y;

	return;
}
// --------------------------------------------------
double plant_2DSOS_RDMD::OneMCTime(double r_d)
{
	double rd;

	rd=RandRand();
	while ((rd <= 0) || (rd >= 1))
		rd=RandRand();

	return -log(rd)/r_d;
}
// --------------------------------------------------
void plant_2DSOS_RDMD::h2mode(){
	int m,n,i,j;
	int index;
	time_t time_begin;
	time_t time_end;
	double dx_square = (pi/LatticeSize)*(pi/LatticeSize);
	std::time(&time_begin);
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D_mode(m,n);
			cov_z1[index] = 0.0;
			cov_z2[index] = 0.0;
			cov_z3[index] = 0.0;
			cov_z4[index] = 0.0;
			for(i=0;i<LatticeSize;i++){
				for(j=0;j<LatticeSize;j++){
					cov_z1[index]+=lattice[i][j]*phi1mn(m,n,i,j);
					cov_z2[index]+=lattice[i][j]*phi2mn(m,n,i,j);
					cov_z3[index]+=lattice[i][j]*phi3mn(m,n,i,j);
					cov_z4[index]+=lattice[i][j]*phi4mn(m,n,i,j);
				}
			}
			cov_z1[index]*= dx_square;
			cov_z1[index] = cov_z1[index]*cov_z1[index];
			cov_z2[index]*= dx_square;
			cov_z2[index] = cov_z2[index]*cov_z2[index];
			cov_z3[index]*= dx_square;
			cov_z3[index] = cov_z3[index]*cov_z3[index];
			cov_z4[index]*= dx_square;
			cov_z4[index] = cov_z4[index]*cov_z4[index];
		}
	}
	std::time(&time_end);
	printf("h2mode() takes %f\n",std::difftime(time_end,time_begin));
}
/************************************************************************/
void plant_2DSOS_RDMD::h2mode_fftw(){
	time_t time_begin;
	time_t time_end;
	std::time(&time_begin);
	
	std::time(&time_end);
	printf("h2mode() takes %f\n",std::difftime(time_end,time_begin));
}
/************************************************************************/
/* Choose and return an initial random seed based on the current time.
   Based on code by Lawrence Kirby <fred@genesis.demon.co.uk>.
   Usage: srand (time_seed ()); */
unsigned time_seed (void)
{
  time_t timeval;	/* Current time. */
  unsigned char *ptr;	/* Type punned pointed into timeval. */
  unsigned seed;	/* Generated seed. */
  size_t i;

  timeval = time (NULL);
  ptr = (unsigned char *) &timeval;

  seed = 0;
  for (i = 0; i < sizeof timeval; i++)
    seed = seed * (0xff + 2u) + ptr[i];

  return seed;
}
/*****************************************************************************/
double ranO(int *idum)
/*Returns a uniform random deviate between 0.0 and 1.0 using rand().*/
/*Set idum to any negative value to initialize or reinitialize the sequence*/
{
	static double y, maxran, v[131];
	double dum;
	static int iff=0;
	int j;

	if (*idum < 0 || iff == 0)
	{
		iff = 1;
		maxran = RAND_MAX+0.0;
		srand((unsigned)time(NULL));
		*idum = 1;
		for (j=0; j<=130; j++)
		{
			dum=rand()+0.0;
		}
		for (j=0; j<=130; j++)
		{
			v[j]=rand()+0.0;
		}
		y=rand()+0.0;
	}
	j= (int)(130.0*y/maxran);
	if (j > 130 || j < 0)
	{
		printf("RanO: This cannot happen.\n");
	}
	y=v[j];
	v[j]=rand();
	return y/(maxran+1.0);
}

/************************************************************************/
double RandRand()
{
	double r_tot, r;
	int i, n;
	static int seed;

	n=3;
	r_tot=0.0;

	for (i=0; i<n; i++)
	{
		r=ranO(&seed);
		if (i > 0)
		{
			r/=pow((RAND_MAX+0.0), (i+0.0));
		}
		r_tot+=r;
	}
	
	return r_tot;
}
double RandN1()
{
	double rdm;

	while((rdm=RandRand())>=1.0) {}

	return rdm;
}
/************************************************************************/
